clear; clc;
import com.comsol.model.*
import com.comsol.model.util.*
model = ModelUtil.create('Model');
model.component.create('comp1', true);
geom1 = model.component('comp1').geom.create('geom1', 3);


geom1.feature.create('blk1', 'Block');
geom1.feature('blk1').set('size', [5, 5, 3]); 
geom1.feature('blk1').set('pos', [0, 0, 0]);  

holeNum = 200;
minHoleRadius = 0.05;       
maxHoleRadius = 0.15; 
cubeSize = [5, 5, 3]; 

for i = 1:holeNum
   
    holeRadius = minHoleRadius + rand() * (maxHoleRadius - minHoleRadius);
    cx = holeRadius + rand() * (cubeSize(1) - holeRadius);
    cy = holeRadius + rand() * (cubeSize(2) - holeRadius);
    cz = holeRadius + rand() * (cubeSize(3) - holeRadius);

    sphereName = ['sph' num2str(i)];
    geom1.feature.create(sphereName, 'Sphere');
    geom1.feature(sphereName).set('r', holeRadius);
    geom1.feature(sphereName).set('pos', [cx, cy, cz]);
    holeVolumes(i) = (4/3) * pi * holeRadius^3;
end

cubeVolume = prod(cubeSize);
totalHoleVolume = sum(holeVolumes);
porosity = totalHoleVolume / cubeVolume;  
fprintf('===================== 空隙率计算结果 =====================\n');
fprintf('立方体总体积：\t\t%.8f\n', cubeVolume);
fprintf('所有孔洞累计体积：\t%.8f\n', totalHoleVolume);
fprintf('多孔介质空隙率：\t%.8f (%.7f%%)\n', porosity, porosity*100);
fprintf('=========================================================\n');

geom1.run();
model.save('3D_porous_medium1.mph');
